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Abstract 

We describe the implementation of analytical Hartree-Fock gradients for pe- 
riodic systems in the code CRYSTAL, emphasizing the technical aspects of 
this task. The code is now capable of calculating analytical derivatives with 
respect to nuclear coordinates for systems periodic in 0, 1, 2 and 3 dimensions 
(i.e. molecules, polymers, slabs and solids). Both closed-shell restricted and 
unrestricted Hartree-Fock gradients have been implemented. A comparison 
with numerical derivatives shows that the forces are highly accurate. 
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I. INTRODUCTION 



The determination of equilibrium structure is one of the most important targets in elec- 
tronic structure calculations. In surface science especially, theoretical calculations of surface 
structures are of high importance to explain and support experimental results. Therefore, 
a fast structural optimization is an important issue in modern electronic structure codes. 
Finding minima in energy surfaces is substantially simplified by the availability of analytical 
gradients. As a rule of thumb, availability of analytical gradients improves the efficiency by 
a factor of order with N being the number of parameters to be optimized. UK's Collab- 
orative Computational Project 3 has therefore supported the implementation of analytical 
gradients in the electronic structure code CRYSTAL&I. This implementation will also be 
valuable for future projects which require analytical gradients as a prerequisite. Another 
advantage of having analytical gradients is that higher derivatives can be obtained with less 
numerical noise (e.g. the 2nd derivative has less numerical noise when only one numerical 
differentiation is necessary). 

CRYSTAL is capable of performing Hartree-Fock and density-functional calculations for 
systems with any periodicity (i.e. molecules, polymers, slabs and solids). The periodicity 
is "cleanly" implemented in the sense that, for example, a slab is considered as an object 
periodic in two dimensions and is not repeated in the third dimension with one slab being 
separated from the others by vacuum layers. The code is based on Gaussian type orbitals 
and the technology is therefore in many parts similar to that of molecular quantum chemistry 
codes. As the density-functional part of the code relies in big parts on the Hartree-Fock 
part, the strategy of the project was to implement Hartree-Fock gradients first. 

The implementation of Hartree-Fock gradients for multicenter basis sets was pioneered by 
Pulayi; the theory had already been derived earlier independentlyil. Meanwhile, analytical 
gradients have been implemented in many molecular codes, and several review articles have 
appeared (see, e.g., references [?|-[T3|). 

Substantial work has also been done in the case of one-dimensional periodicity: Hartree- 
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Fock gradients with respect to nuclear coordinates and with respect to the lattice vector have 
already been implemented in codes periodic in one dimensionlli^El. Moreover, correlated 
calculations based on the MP2 scheme0'0 and MP2 gradients^ have been coded. Also, 
density functional gradients have been implemented^^. Even second derivatives at the 
Hartree-Fock level have meanwhile been coded0. 

The implementation of Hartree-Fock gradients with respect to nuclear coordinates in 
CRYSTAL is to the best of our knowledge the first implementation for the case of 2- and 
3-dimensional periodicity. The aim of this article is to describe the implementation of the 
gradients in the code, with an emphasis on the technical aspects. Therefore, the article is 
supposed to complement our first article on the purely theoretical aspectsS. An attempt of 
a detailed description is made; however, as the whole code is undergoing constant changes, it 
can not be too detailed. For example, it did not seem advisable to give any variable names 
because they have already undergone major changes after the code moved to Fortran 90 
with the possibility of longer variable names. 

The article is structured as follows: In section we give a brief introduction to Gaussian 
and Hermite Gaussian type basis functions. The definition of the density matrix is given in 
section |TT[ The individual integrals, their derivatives, and details of the implementation are 
discussed in section Formulas for total energy and gradient are given in section 0. The 



structure of the gradient code is explained in section VI, followed by examples in section 



VII and the conclusion. 



II. BASIS FUNCTIONS 

Two sets of basis functions are relevant for CRYSTAL: firstly, unnormalized spherical 
Gaussian type functions, in a polar coordinate system characterized by the set of variables 
(|r|, i), If), and centered at A. They are defined as 

S{a,r — A,n,l,m) = \f — A\ ^P|'"'(cos ^d) exp^imip) exp{—a\f — A\ ) (1) 
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with Pp' being the associated Legendre function. CRYSTAL uses real spherical Gaussian 
type functions defined as 

R{a, r- A, n, I, 0) = S{a, f-A,n, 1, 0) 
R{a,f— A, n, I, \m\) — Re S{a,f— A, n, I, \m\) 
R{a,f— A,n,l,—\'m\) = Im S{a, f — A,n,l,\m\) 

This is in the following denoted as (j)fj,{a, f— A^j^, n, I, m) — Nf^R{a, r — A^, n, I, m), with 
the normalization A*"^. fj, is an index enumerating the basis functions in the reference cell 
(e.g. the primitive unit cell). In fact, CRYSTAL uses only basis functions with quantum 
number n — and angular momentum 1—0,1 or 2 (i.e. s, p or d functions). 

The exponents are defined by the user of the code. A huge amount of basis sets for 
molecular calculations is available in the literature and on the world wide web; also for 
periodic systems a large number of basis sets has been published. Molecular basis sets can, 
with a little effort, be adopted for solid state calculations. High exponents which are used to 
describe core electrons do not have to be adjusted, but exponents with low values (e.g. less 
than 1 Oq ^, with gq being the Bohr radius) should be reoptimized for the corresponding solid. 
Very diffuse exponents should be omitted because they cause linear dependence problems 
in periodic systems. 

A second type of basis functions, which CRYSTAL uses internally to evaluate the inte- 
grals, is the Hermite Gaussian type function (HGTF) which is defined as: 

A(,,r-A.t.u.v)={^)\^)\^)\M-l\r-m (2) 

CRYSTAL uses the McMurchie-Davidson algorithm to evaluate the integrals. The basic 
idea of this algorithm is to map the product of two spherical Gaussian type functions on 
two centers onto a set of Hermite Gaussian type functions at one center. 

S{a,f— B, h, I, rn)S{a,f — A, n, I, m) — 

E{n, I, m, n, I, m, t, u, v)A{j, f— P, t, u, v) (3) 

t,u,v 
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with 7 = a + a and P = 

The starting point £;(0, 0, 0, 0, 0, 0, 0, 0, 0) = exp(-^|fi - Ap) is derived from the 
Gaussian product rule@: 

expf— a r — A\ ) exp[—a\r — B\ ) = exp -\n — A\ exp — [a + a) r z — 

\ a + a J \ a + a 



As indicated in section all the integrals can be expressed with the help of the 
coefficients E(n,l,m,n,l,m,t,u,v)^^^. These coefficients are generated by recursion 
relationsSii. They are zero for the case t + u + v > 2n + 2h + I + I and for all nega- 
tive values of t, u or v. CRYSTAL uses only basis functions with n = 0. Therefore, there are 
(;+;+i)(;+^+2)(;+;+3) coefficients E(0,l,7fi,0,l,m,t,u,v) for fixed values of l,m,l,m. As the 
maximum angular quantum number is / = 2, this results in 25 possible combinations of m 
and rh. Therefore, the maximum number of coefficients is 25 x 35 = 875. These coefficients 
are pre-programmed in the subroutine DFAC3. Pre-programming is the fastest possible way 
of evaluating these coefficients which is important because this is one of the key issues of the 
integral calculation. On the other hand, the code has become inflexible as no iJ-coefficients 
are available for higher quantum numbers. 

Derivatives of Gaussian type functions are again Gaussian type functions. Therefore, 
the evaluation of gradients is closely related to the evaluation of integrals. In a similar way 
as all the integrals can be expressed with the help of coefficients E, all the derivatives of 
the integrals can be expressed with the help of coefficients for the gradients, G';^,G'^,G^. 
These G-coefficients can be obtained with recursion relations derived by Saunder^B. The 
recursions are similar to the ones for the ^^-coefficients. However, as the existing subroutine 
DFAC3 cannot compute the G-coefficients, the recursions were newly coded. This has in 
addition the advantage that, by small modifications of the new subroutines, E'-coefficients 
for higher quantum numbers than 1 = 1 = 2 can now be computed by recursion. There are 
three sets of G-coefficients because of the three spatial directions. The G-coefficients are 
zero for the case t + u + v > 2?7, + 2fi + / + / + l and for all negative values of t, u or v. 
This means that for a maximum quantum number of / = 2, there are 3x5x5x 56 = 4200 



coefficients. Three other sets of G-coefficients are necessary because of the second center. 
However, the sets on the second center are closely related to the sets on the first center and 
can be derived from them in an efficient wayi'ii. 

III. DENSITY MATRIX 

After solving the Hartree-Fock equationsil, the crystalline orbitals are linear combina- 
tions of Bloch functions 

^i{f,k) = J2af,^{k)^^J^{f,k) (5) 
which are expanded in terms of real spherical Gaussian type functions 

^^(f, k) = N^Y. ^(«> r-A,- g, n, I, my'^ (6) 

9 

The sum over g is over all direct lattice vectors. 

In the case of closed shell, spin-restricted Hartree-Fock, the spin-free density matrix in 
reciprocal space is defined as 

P,,{k) = 2j2a^,{k)aU{k)e{eF - e^)) (7) 

i 

with the Fermi energy ep and the Heaviside function B; z is an index enumerating the 
eigenvalues. 

In the case of unrestricted Hartree-Fock we use the notation 

^]{r,k)=Y.aUk),lj,{r,k) (8) 

and 

m\{r,k) = Y.aim^,{r,k) (9) 

for the crystalline orbitals with up and down spin, respectively. We define the density 
matrices 
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PUk) = Y.'^UkKUk)Qiep - e]{k)) (10) 

i 

for up spin and 

Pluik) = Y.ciUk)aMe{eF - e\{k)) (11) 

% 

for down spin. In the following, P^,^ refers to the sum Pj^^ + P"^^ in the UHF case. 
The density matrices in real space -P^Oi^^', P^iOiycf P^oug obtained by Fourier transfor- 
mation. 



IV. INTEGRALS AND THEIR DERIVATIVES 

The calculation of the integrals is fundamental to all quantum chemistry programs. 
CRYSTAL uses two integral packages: a package derived from GAUSSIAN70il is the default 
for calculations when only s and sp shells are used; alternatively Saunders' ATMOL Gaussian 
integral package can be used and it must be used for cases when poi d functions are involved. 
The implementation of gradients has been done with routines based on the ATMOL package. 
This is not a restriction, and it is possible to use routines based on GAUSSIAN70 for the 
integrals and routines based on ATMOL for the gradients. 

The calculation of the integrals is essentially controlled from MONMAD and MONIRR 
for one-electron integrals and from SHELLC or SHELLX for the bielectronic integrals. 
SHELLC is used in the case of non-direct SCF, i.e. when the integrals are written to 
disk and read in each cycle. SHELLX is the direct version when the integrals are computed 
in each cycle without storing them on disk. The direct mode is the preferred one when the 
integral file is too big or when input /output to disk is too slow. The gradients are com- 
puted only once after the last iteration, when convergence is achieved. Therefore, a direct 
implementation of gradients has been done. 

One of the bottlenecks of the CRYSTAL code is the restriction to a highest quantum 
number of / = 2, i.e. the code can only cope with s, p, sp and d functions, but not with 
basis functions with higher angular momentum. Introducing gradients, however, is similar 
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to increasing the quantum number from d to f for the corresponding basis function. This 
means that many subroutines had to be extended to higher quantum numbers, and array 
dimensions in the whole code had to be adjusted. 



A. One-electron integrals 



In this section we summarize the appearing types of integrals and the corresponding gra- 
dients. We restrict the description to the x-component of the gradient; y- and z-component 
can be obtained in similar way. Note that the integrals depend on the dimension because 
of the Ewald scheme used. Therefore, there are four different routines for the one-electron 
integrals for the case of 0,1,2 and 3-dimensional periodicity: CJATO, CJATl, CJAT2 and 
CJAT3. Similarly, four gradient routines have been developed which have been given the 
preliminary names CJATOG, CJATIG, CJAT2G and CJAT3G. These routines calculate all 
the one-electron integrals except for the multipolar integrals which are computed in POLIPA 
(with the corresponding gradient routine OLIPAG). 



1. Overlap integral 



The basic integral is the overlap integral: 




(12) 
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E{n, I, m, n, I, m, 0, 0, 0) 




The x-component of the gradient with respect to center is obtained as 



d 



(13) 



OA, 
d 



dA, 



dA, 
d 
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t,u,v 



G^»(..[,*.M,m.O.O.O)l^| (14) 



Equation IT3 thus 



defines the coefficients G^"; similarly the coefficients Gy^',G^^',G^'',Gy'',G^" can be de- 
fined. 

In the following, we use the identity 

2. Kinetic energy integrals 



In equation 15, the expression for the kinetic energy integrals for the case of spherical 



Gaussian type functions is reiteratedE^: 



T - ^ 



(f)^{a,f- A^,n, l,m) (-^Af^ (j)y{a,r - - g,nj,m)dh = 
«(4n + 2( + 3) / ^,(a, f - i„ n, ( , m)Ma, ? - X - 9. «, /, m)d'r - 
—n{2n + 2/ + 1) j ^ E{n, l,ih,n — 1, /, m, t, u, f )A(7, r — P, t, u, v)(i^r + 

t,u,v 

a{An + 2/ + 3) J ^ E{n, [, m, n, /, m, t, f )A(7, r — F, t, -u, v)d^r — 

t,u,v 

2a^y ^ i?(fi, [, m, n + 1, /, m, t, M, t;)A(7, r — P, t, M, t')d'^r (15) 



t,u,v 



The x-component of the gradient is therefore: 
d 



(2n + 2/ + l) y ^ G'^''(n,[,m,n - l,/,m,t,M,t;)A(7,f- P,t,M,w)d^r + 



-n 



t.U'U 
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a{An + 2/ + 3) / ^ G^^in, /, m, n, /, m, t, m, 'y)A(7, f - P, t, u, v)dh - 

•' t,u,v 

2o? J G^^in, l,m,n + 1, /, m, t, u, v)A{'y, r- P, t, u, v)dh (16) 



t,u,v 



As CRYSTAL uses spherical Gaussian type functions with n = 0, this reduces to 



dA^,, ^^^^ 

3 
2 

a(4n + 2/ + 3)G'^'*(0, /, m, 0, /, m, 0, 0, 0) - 



2 j a^G^''(0, /, m, 1, /, m, 0, 0, 0) (17) 

Exphcit differentiation with respect to the other center Ai, is more difficuh because the 
kinetic energy operator apphes to that center. However, the differentiation can easily be 
avoided by applying translational invariance: 



3. Nuclear attraction integrals 
The nuclear attraction integrals are defined as 

^t^oug = - H^a I 4>i,{ot, r- A^, fi, I, 'm)A{f- Aa)(pi.{a, f- A^-g,n, I, m)dh = 

- II ^« / J2 ^(^' ^~ ^ ^' ^' r- P, t, u, v)A{f- Aa)dh (19) 

a •' t,u,v 

where A is the Coulomb potential function in the molecular case, the Euler-MacLaurin 
potential function for systems periodic in one dimension!!. Parry's potential functionll for 
systems periodic in two dimensions, and the Ewald potential function for systems periodic 
in three dimensionsSSB. The summation with respect to a runs over all nuclei of the 
primitive unit cell. 

The x-component of the partial derivative with respect to the center A^^^ is obtained as: 
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-E^"/ Y.'^x''^nJ,m,n,l,m,t,u,v)K{-i,r- P,t,u,v)A{r- Aa)dh (20) 

a t,u,v 

In the same way, the partial derivative with respect to A^^^ is obtained. The partial 

derivative with respect to the set of third centers Aa is obtained by translational invariance: 

for each center Aa, there is a derivative with value 
a d_ 

4. Multipolar integrals 
The electronic charge density is expressed with a lattice basis as: 

P(^) = - J2 Kgt.d'PAo!, f- Af,, h, I, m)(j)y{a,r- A^ - g, n, I, m) (21) 
Then, the Ewald potential due to this charge density is given by: 

^"^""{p-r) = j A{f-f^)p{f^)Sr' (22) 

The Ewald energy of the electons (i.e. the Ewald energy of the electrons in the primitive 
unit cell with all the electrons) is obtained as 

E=^l I p{r)A{f- f^)p{f^)Srd^r' (23) 

For efficiency reasons, the calculation of the Ewald potential is done approximatively. 
A multipolar expansion up to an order L is performed for the charge distribution in the 
long range. Therefore, the electrons do not feel the Ewald potential created by the correct 
charge distribution, but the Ewald potential created by the multipolar moments. It is thus 
necessary to compute the multipolar moments of the charge distribution which are defined 
as 

riT[Pc, A,) = J Pc(rlXr(r - A,)d'T (24) 

with XJ^ being regular solid harmonic and the charge Pc{r) defined as 
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Pcir} = - Kg^^6'Pt^i'^^'^-^t^^'^^'^^^)M(^^r-A^-g,nJ,m) = 

g,^&c,v t,u,v 

c is an index for the shell. The total electronic charge p{r) is thus obtained by summing 
over all shells c: 

P(0=EMr1 (26) 

c 

In CRYSTAL, the multipole is located at center ^4^ and therefore it is convenient to take 
the derivative with respect to center A^, and obtain the derivative with respect to A^^ by 
translational invariance. The expression computed for the gradients is thus 

~ J2 ^ugt^o J ^- ^1^^ I m)(f)^{a, f- A^ - g, n, /, m)X'^{r- A^)) d^r = 

g,fi&c,v t,u,v 



5. Field integrals 

If the electronic charge distribution is approximated with an expansion up to the maxi- 
mum quantum number L, the Ewald potential of this model charge distribution is obtained 

as 

^e^^pmodel. ^ ^ r") = ^ E E VriPc, A)Zr{t)A{r- A) (28) 

c c 1=0 m=-l 

with Zl^{Ac) being the spherical gradient operator in a renormalized form@. The model 
charge distribution is expressed as 

pr'^'if^ = E E C(Pc; AW A r-) (29) 

1=0 m=-l 

and 

(5r(i'c,r) = lim Zr(ie)A(a,r-A„ 0,0,0) (30) 



12 



The integral of the electronic charge distribution and the Ewald potential function is 
required which gives rise to the field integrals which are defined as follows: 



Zl^{Ac) / E{n, I, fh, n, I, m, t, u, f )A(7, r — P, t, u, v) 



t,u,v 



A{f-A, 



Air -A,) 



pen 

E 

n 
pen 

E- 



Ir — Ac — n\ 



r- Ac 



n\ 



d^r = 
d^'r (31) 



The term 



A{r-Ac)-Y. 



pen 



instead of A{f—A()j appears because the multipolar 



\r-Ac-n\ 

approximation is only done for the charge distribution in the long range. The penetration 
depth pen is a certain threshold for which the integrals are evaluated exactly@^. 

For the gradients, the derivative with respect to all the centers is needed. The partial 
derivative with respect to is obtained as 



Ml 



Zr{Ac) / ^ G^(n,/,m,n,/,m,t,w,v)A(7,f-F,t,w,^;) y4(r-^)-^ 



t,u,v 



pen 



,3 \r — Ar — n\ 



d^r (32) 



In similar way, the partial derivative with respect to center A^, is computed. Finally, the 
partial derivatives with respect to the centers A^. are obtained from translational invariance. 



6. Spheropole 

This term arises because the charge distribution is approximated by a model charge 
distribution in the long rangeil: 

The calculation of the Coulomb potential <|)COM«('p^_pmo(ie«. ^ restricted to contributions 
from those charges inside the penetration depth pen. The use of the Coulomb potential 
^coui (^p^_ puiodei r) instead of the Ewald potential $'''"(pc - p™"'^'''; f) is correct, if pc - p™°^''' 
is of zero charge, dipole, quadrupole and spherical second momentil. However, this condition 
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leads to a correction in the three-dimensional caseBElil: although the difference pc — p^"'^^^ 
has zero charge, dipole and quadrupole moment, it has in general a non-zero spherical second 
moment Qc- Therefore, the potential must be shifted by Q defined as: 

Q = Ea = E^/ - pr'^'(r-))|rTd^r (34) 



Three types of contributions are obtainedE^: zero, first and second order HGTFs. They 
have to be combined with the corresponding i?-coefficient. For the zeroth order, a contri- 
bution of 

E{h, I m, n, I, m, 0, 0, 0) (^^ + [A^ - P) 
is computed. The derivative is therefore 

^ ' E{fiJ,m,nJ, m,0, 0,0) {^ + {A^-P)A] (35) 



E{h, I, m, n, I, m, t, u, v)A{'j, f — P, t, u, v) 



To obtain the derivative — -E(n, /, m, n, I, m, 0, 0, 0), we use the identity 
d 



'-'^fl,X \t,u,v 

E{h, /, m, n, I, m, t, u, f )— A(7, f — P,t + l,u,v) + 

t,u,v T 

d ^ ~ ^ 

E A(7, f - P, t, w, v)— — E{n, I, m, n, I, m, t, u, v) = 



t.u.v 



G^^in, /, m, 77,, /, m, t, u, f )A(7, r — P, t, u, v) (36) 

which gives 

E{n, I, m, n, I, m, t, u, v) = G^''{n, I, m, n, I, m, t, u, v) — —E{h, I, m, n, l,m,t — 1, u, v) (37) 



oA^^x 7 

A similar operation is necessary for the components with i?(n, /, m, n, m, 1, 0, 0), 
£'(fi, /, m, r?,, /, m, 0, 1, 0) and i?(n, /, m, tt,, /, m, 0, 0, 1) (first order HGTFs) which are mul- 
tiplied with prefactors 2(F^. — A^^x)^ 2(Py — A^^y) and 2(Pj — A^,^), respectively. Fi- 
nally, derivatives of the products of E{n, I, 'm,n, I, m, 2, 0,0), E{n, I, in, n, I, m, 0,2,0) and 
E{n, I, fh, n, I, m, 0, 0, 2) (second order HGTFs) with 2 are required. 
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B. Bielectronic integrals 

We define a bielectronic integral as 

^Ougrnan+h 

f-A^, ni, /i, mi)0^(a2, r-A-g, n2, k, ^2) 

(f)r{c(3,f'—AT- — n, 72-3, /s, m3)0o-(a4, r'—Aa- — n — h, 714, ^4, m4)d^r d'^r' = 

E{ni,li,mi,n2,l2,m2,t,u,v) ^ E(n3, ^3, ma, n^, k, m^, t' , u', v') [t, u, v\ ^ u' , v'] (38) 

t,u,v t'u'v' 



The expression [t, u, v\ ,^}^n \t', u', v'] is defined a^iH 



\r—r'\ 

[t, f I p ^ u', v'] = 

A(7, f - P, t, u, i;)^r^A(7 ',f' -P ', t,' u', v')dh ' (39) 

The partial derivative with respect to A^^x is obtained as 
d 



-B 



dA pLOvgrnan+h 

^ G^''{ni,li,mi,n2,l2,m2,t,u,v) ^ E{n3,l3,m3,n4,,U,m4,t' ,u ,v')[t,u,v\ ^ ^ ^ \t',u',v'] (40) 

t,u,v t',u',v' I'" ^ I 

Similarly, gradients with respect to the other centers are obtained. One of the gradients 
can be obtained by translational invariance if the other three gradients have been computed. 

In the context of periodic systems, it is necessary to perform summations over the lattice 
vectors g, h, n. We define a Coulomb integral as follows 

pen 

(41) 

n 

Similarly, we define an exchange integral as follows: 

fiOrnvgcrn+h 

(42) 
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V. TOTAL ENERGY AND GRADIENT 



A. Total energy 

The correct summation of the Coulomb energy is the most severe problem of the total 
energy calculation. The individual contributions to the Coulomb energy, such as for example 
the nuclear-nuclear interaction, are divergent for periodic systems. Thus, a scheme based 
on the Ewald method is used to sum the individual contributions^. The total energy is 
then expressed as the sum of kinetic energy S'^^", the Ewald energies of the nuclear-nuclear 
repulsion E^^ , nuclear-electron attraction £'coM«-n«c^ electron-electron repulsion E'^°'^^~^\ 
and finally the exchange energy E^^'^^~^K 

^total ^kinetic _j_ ^NN _j_ ^coul— nuc _j_ ^coul— cl _j_ ^oxch— cl 

— P - -^ F^^ 

g,^^,u T o- c l=Q m=-l 

1 1 

— - ^ p^^ X- — - T^p^^y^ p^^ x ^ - - (43) 

2 ^ ugfj.6 ^ ahrO l-tOugrOah 2 ^ i^gf^O ^ crhrO t-tOugrOcrh \ ) 



B. Gradient of the total energy 

The force with respect to the position of the nuclei can be calculated similarly to the 
molecular caseiB. The derivatives of all the integrals are necessary, and the derivative of the 
density matrix is expressed with the help of the energy- weighted density matrix. The full 
force is obtained as: 

^^total 



dAi 



ugfiO 



g,fi,v dAi 



16 



1 



d_ 



S.. 



2-n 



V("2, ^- A^, n2, h, m2)(t)u{oii,r - Ay - g,ni,h, mi)A{f- Aa)dh 

h,a,T£c 



d 



dA 



- 0r(a2, r - yl^, n2, 12, m2)(f)aiai, f- A^- h, ni, /i, mi) 

+ E E / 0r(«2,r'-l,,n2,/2,m2)(/).(ai,f'-X-K,ni,/i,mi)Xr(f'-le)dV5™(l^ 
^ dC. 



hrO 



T,<T dAi 

L I 



+ E E E E ^ahrQ 

C l=Q m=-l h^^izc^^ 

d 



dA, 
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(44) 



The last addend is the energy weighted density matrix; the integral is over the first 
Brillouin zone. 



VI. STRUCTURE OF THE GRADIENT CODE 

The present structure of the gradient code is indicated in figure |l]. The first step is to 
compute the gradient of the Ewald energy of the nuclei in subroutine GRAMAD (the Ewald 
energy is computed in ENEMAD). The control module TOTGRA then first calls routines to 
compute the gradient of the bielectronic integrals (labeled with SHELLXV as these routines 
will change their structure). The subroutine SHELLXV calls subroutines which explicitly 
compute the derivatives of Coulomb and exchange integrals, and multiplies the gradients a 
first time with the density matrix. Back in TOTGRA again, the second multiplication with 
the density matrix is performed. The next step is to compute the derivatives of the multi- 
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poles (MONIRG) and to compute the energy weighted density matrix (PDIGEW). Then, 
the gradients of the one-electron integrals are computed (CJATOG, CJATIG, CJAT2G or 
CJAT3G, depending on the dimension). The field integrals and their gradients are now 
multiplied with the multipolar integrals and their gradients, and a multiplication with the 
density matrix is performed. This concludes the calculation of the gradients. 

The structure has been simplified to focus on the most important parts. In addition, as 
already mentioned, the code will undergo changes during the optimization process so that 
a too detailed description seems to be unadvised. 
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FIGURES 

FIG. 1. The present structure of the gradient code. The left column describes the purpose 
of the routines, the middle column gives the names of the corresponding routines, and the right 
column gives the name of the routines in the energy code. One arrow indicates that the routine is 
a subroutine, two arrows indicate that it is a subroutine called from a subroutine. 
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VII. EXAMPLES 



In tables |, |T| and |T|, we give examples of the accuracy of the gradients. First, in table 
|, a chain of NiO molecules is considered, with ferromagnetic ordering (all the Ni spins up) 
and with antiferromagnetic ordering (nearest Ni spins are antiparallel). The oxygen atoms 
are moved by 0.01 A from their equilibrium positions which results in a non-vanishing 
force. The agreement between numerical and analytical gradient is better than 0.0001 
As we discussed in our first articled, the agreement can be improved by using stricter 
"ITOL" -parameters (these are parameters which control the accuracy of the evaluation of 
the integrals!). Indeed, when increasing these parameters, the agreement further improves 
up to an error of less than 10~^ — . 

In table a LiF layer with a lattice constant of 5 A is considered with one atom being 
displaced from its equilibrium position. The forces agree to 2 x 10^^-^ when default ITOL 
parameters (6, 6, 6, 6, 12) are used. 

Finally, in table |T|, a three-dimensional, ferromagnetically polarized NiO solid is con- 
sidered. When displacing the oxygen ions, the forces agree to better than 2 x 10^^-^. 

As a whole, the accuracy is certainly very high and can further be improved by applying 
stricter cutoff (ITOL) parameters. 



VIII. CONCLUSION 

In this article, we described the implementation of analytical gradients in the code CRYS- 
TAL. In its present form, the code is capable of computing highly accurate Hartree-Fock 
gradients for systems with 0,1,2 and 3-dimensional periodicity. Both closed-shell restricted 
Hartree-Fock as well as unrestricted Hartree-Fock calculations can be performed. 

A first step of improving the efficiency of the code has been completed with the coding of 
gradients for the bipolar expansion, and a further enhancement of the efficiency will be one 
of the future directions. Of highest importance is the implementation of symmetry which 
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will lead to high saving factorsED. Other targets are the implementation of gradients with 
respect to the lattice vector, an extension to metallic systemsil, and the implementation of 
density functional gradients. 
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TABLES 

TABLE I. Ferromagnetic (FM) and antiferromagnetic (AFM) NiO chain (i.e. a chain with 
alternating nickel and oxygen atoms). The distance between two oxygen atoms is chosen as 5 A. 
The force is computed numerically and analytically with the oxygen atoms being displaced. A 
[5s4p2c(| basis set was used for nickel, and a [4s3p] basis set for oxygen. 



magnetic 
ordering 

FM 
FM 
AFM 

AFM 



ITOL 
parameter 

6 6 6 6 12 
8 8 8 8 14 
6 6 6 6 12 

8 8 8 8 14 



displacement 
of oxygen 
in A 
0.01 
0.01 
0.01 
0.01 



analytical derivative 
(x-component) 
Eh/ao 
0.001274 
0.001246 
0.001276 
0.001250 



numerical derivative 
(x-component) 
Eh/ao 
0.001188 
0.001249 
0.001191 
0.001252 



TABLE II. Forces on the atoms of a LiF layer when one of the atoms is displaced from its 
equilibrium position. A [AsSp] basis set was used for the fluorine atom and a [2slp] basis set for 
the lithium atom. Default ITOL parameters were used. 



atom 

F at (0.5 A, A, A) 
Li at (2.5 A, A, A) 
F at (2.5 A, 2.5 A, A) 
Li at (0 A, 2.5 A, A) 



analytical derivative 
(x-component) 
0.001379 
-0.020731 
0.010384 
0.008969 



numerical derivative 
(x-component) 
0.001400 
-0.020726 
0.010376 
0.008950 
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TABLE III. Ferromagnetic NiO in an fee structure at a lattice constant of 4.2654 A. We com- 
pare numerical and analytical derivatives when moving the oxygen ion parallel to the x-direction. 
Default ITOL parameters were used, the basis sets are the same as in table |. 



displacement of oxygen 

in A 
0.01 
0.02 
0.03 
0.04 
0.05 



analytical derivative 
(x-component) 
Eh/ao 
0.001499 
0.002939 
0.004387 
0.005857 
0.007352 



numerical derivative 
(x-component) 
Eh/ao 
0.001485 
0.002925 
0.004378 
0.005847 
0.007346 
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